.  REPORT  DOCUMENTATION  PAGE _  QMBZAomeoi8B 

PubiiQ^epcrting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducinq 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway  Suite  1204  Arlington  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS.  y 

1,  REPORT  DATE  (DD-MM-YYYY)  [~2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  To) - 

25-06-2003 _ Technical  Paper _ 

4,  TITLE  AND  SUBTITLE  Ra  nnNTRArvr  miimrpo 


Technical  Paper 


3.  DATES  COVERED  (From  -  To) 


Flexible  Three-Dimensional  Modeling  of  Electric  Thrusters  in  Vacuum  Chambers 


5a.  CONTRACT  NUMBER 
F04611-99-C-0025 
5b.  GRANT  NUMBER 


5c,  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S)  ~  - - 

M.R.  Gibbons,  D.E.  Kirtley  (ERC);  D.B.  VanGilder,  J.M.  Fife  (AFRL/PRSS) 


5d.  PROJECT  NUMBER 
4S47 

5e,  TASK  NUMBER 
0052 

51.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


ERC 

10  East  Saturn  Blvd. 

Edwards  AFB  CA  93524-7680 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/PRS 
5  Pollux  Drive 

Edwards  AFB  CA  93524-7048 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S 
NUMBER(S) 

AFRL-PR-ED-TP- 2003-173 


Approved  for  public  release;  distribution  unlimited. 

13.  SUPPLEMENTARY  NOTES  ~ 

For  presentation  at  the  AIAA  Joint  Propulsion  Conference  in  Huntsville,  AL,  20-23  July  2003 

14.  ABSTRACT  '  * - 1 - 


15.  SUBJECT  TERMS 


20030812  157 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 


b.  ABSTRACT  I  c.  THIS  PAGE 


Unclassified  Unclassified  Unclassified 


17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Leilani  Richardson 

A 

9 

19b.  TELEPHONE  NUMBER  (include  arm 
code) 

(661)275-5015 

Prescribed  by  ANSI  Std.  239,18 


Flexible  Three-Dimensional  Modeling  of  Electric  Thrusters  in 

Vacuum  Chambers 


M.R.  Gibbons  and  D.  E.  Kirtley 
ERC,  Inc. 

Edwards  AFB,  CA 

D,  B.  VanGilder  and  J,  M.  Fife 
U.S.  Air  Force  Research  Laboratory 
Edwards  AFB,  CA 


A  national  team  of  researchers  is  developing  a  software  system,  COLISEUM,  which  performs  calculations  of 
plasma  propagation  and  interaction  with  arbitrary  3-D  surfaces.  The  applications  of  COLISEUM  are  wide- 
ranging,  but  include  simulating  engine  test  configurations  inside  vacuum  chambers  and  predicting  sputtering 
and  re-deposition  on  spacecraft  surfaces.  COLISEUM  allows  users  to  easily  define  complicated  3-D  geometries 
using  off  the  shelf  CAD  software,  then  select  from  a  set  of  plasma  expansion  models  of  varying  fidelities  and 
numerical  complexity  to  perform  the  solution.  Once  the  object  surfaces  are  created,  the  user  can  run  different 
types  of  simulations  for  the  same  geometry.  With  this  system,  low  fidelity  models  can  be  used  to  verify  the 
geometry  and  boundary  conditions,  and  to  obtain  first-order  predictions.  Higher  fidelity  models  are  then  used  to 
obtain  more  accurate  predictions  with  greater  cost  in  computation  time.  This  paper  describes  the  simulation  of  a 
Hall  thruster  firing  inside  a  vacuum  chamber,  and  the  execution  of  two  types  of  simulations  for  the  same 
geometry;  one  in  which  the  equilibrium  sputter  and  deposition  rates  are  calculated,  and  another  in  which  the 
Hall  thruster  plume  expansion  is  calculated.  The  sputter  and  deposition  calculations  are  made  with  an  iterative 
ray  tracing  algorithm  which  takes  into  account  re-sputtering  of  deposited  material.  The  plume  expansion 
calculations  are  performed  with  a  particle-in-eell  (PIC)  algorithm  which  includes  wall  collisions  and  wall 
recombination.  Techniques  are  incorporated  that  decrease  the  computation  time  required  for  the  PIC 
simulation  to  relax  to  steady  state. 


Introduction 

High-energy  HET  exhaust  ions  may  erode  (sputter)  surfaces 
on  which  they  impinge.  In  addition,  this  sputtered  material 
may  be  re-deposited  on  other  spacecraft  surfaces.  These 
issues,  and  others,  such  as  electromagnetic  interference  and 
spacecraft  charging,  cause  some  concern  for  spacecraft 
designers  who  want  the  maneuverability  EP  offers  but  do  not 
want  increased  risk. 

Efforts  are  underway  to  accurately  quantify  some  of  the  risks 
associated  with  integration  of  EP  with  spacecraft,  including 
surface  erosion  and  re-deposition.  Work  has  been  done  to 
computationally  model  expansion  of  HET  plumes.1 
Additionally,  Gardner  et  al.  have  developed  Environment 
Work  Bench  (BWB),  a  program  that  calculates  sputtering  of 
spacecraft  surfaces  by  superimposing  pre-computed  EP 
plumes  onto  spacecraft  geometries,2^  However,  existing 
programs  do  not  self-consistently  calculate  the  plume 
expansion  with  3-D  surface  sputtering  in  a  usable,  flexible 
way. 

The  Air  Force  Research  Laboratory  is  leading  development  of 
a  new  software  package  named  COLISEUM,  which  is 
capable  of  self-consistently  modeling  plasma  propagation  and 
interactions  with  arbitrary  3-D  surfaces.  Three  important 


requirements  have  been  placed  on  COLISEUM:  It  must  be 
USABLE,  FLEXIBLE,  and  EXPANDABLE. 

USABLE  means  a  typical  engineer  is  able  to  set  up  and  ran  a 
typical  low-fidelity  case  in  less  than  one  day  with  less  than 
three  days  training. 

FLEXIBLE  means  COLISEUM  is  able  to  simulate  at  least 
three  important  cases:  a)  a  single  spacecraft,  b)  multiple 
spacecraft  in  formation,  and  c)  laboratory  conditions  (e.g.  the 
interior  of  a  vacuum  test  facility).  Simulating  laboratory 
conditions  is  very  important  for  two  reasons.  First,  since 
there  is  very  little  on-orbit  data  for  EP  thrusters,  ground-based 
tests  must  be  relied  upon  for  the  bulk  of  code  validation. 
Second,  by  modeling  the  laboratory  conditions,  COLISEUM 
can  help  engineers  interpret  lab  measurements. 

In  addition  to  being  able  to  simulate  multiple  geometries, 
COLISEUM  is  flexible  in  its  use  of  plasma  simulation 
algorithms.  Problem  set-up  and  geometry  definition  is 
preformed  once.  Then,  the  user  may  select  from  a  set  of 
interchangeable  plasma  simulation  algorithms  to  perform  the 
solution.  If  fast  execution  is  desired,  a  low-fidelity  technique 
can  be  selected  such  as  ray  tracing.  For  higher  fidelity  (at  the 
cost  of  longer  run-time),  something  like  Particle-In-Cell  (PIC) 
can  be  used. 
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*  Prescribed  plume  (AFRL) 

*  Ray  tracing  (AFRL) 
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.  .  *  Colloid  Engine  Sim  (Advatech) 

Simulation  .  PIC-lmrnersed-Hybrid-DSMC{¥T) 
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*  PIC-PIC-DSMC 
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Output 
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Fig,  1,  Architecture  for  using  various  interchangeable  plasma  simulation  techniques  with  the  same  3- 
D  surface  geometry. 


EXPANDABLE  means  COLISEUM  can  be  easily  expanded 
to  incorporate  new  plasma  simulation  algorithms,  new 
capabilities,  or  improved  efficiency.  Furthermore,  as  new 
plasma  simulation  algorithms  are  added,  old  ones  will 
continue  to  function. 

Code  Architecture 

Fig.  1  shows  how  the  COLISEUM  Application  Programming 
Interface  (API)  works  with  a  set  of  various  interchangeable 
plasma  simulations  (applications).  In  general,  the 
COLISEUM  API  can  be  viewed  as  a  framework  in  which  3-D 
plasma  simulations  can  be  quickly  integrated.  Common 
calculations  (such  as  those  related  to  surfaces,  material 
properties,  and  flux  sources)  are  standardized,  grouped,  and 
provided  as  a  resource  library  (data  and  subroutines)  to  each 
simulation.  This  resource  library  takes  the  form  of  a  .lib 
file  that  users  link  with  their  set  of  plasma  simulation 
routines. 

Plasma  simulation  modules  are  the  primary  functional 
components  of  COLISEUM.  They  calculate  plasma 
propagation  of  matter  on  the  volume  domain.  They  contain 
algorithms,  such  as  ray  tracing,  fluid,  PIC,  DSMC,  or  hybrids 
thereof,  which  perform  a  solution  subject  to  pre-set  boundary 
conditions.  Plasma  simulation  modules  are  interchangeable. 
They  all  conform  to  the  COLISEUM  API  which  is  formalized 
in  the  Interface  Control  Document  (ICD). 

The  COLISEUM  resource  library  functions  support  tasks 
common  to  all  types  of  plasma  simulations.  They  handle 
boundary  conditions,  and  provide  support  to  plasma 
simulation  modules. 

The  purpose  of  this  modular  design  is  to  give  COLISEUM 
flexibility  and  expandability.  A  large  number  of  plasma 
simulation  modules  are  desired  to  allow  flexibility  in  solving 
a  variety  of  different  problems.  The  ICD  Is,  therefore,  very 
important,  because  it  describes  for  authors  of  plasma 


simulation  modules  a)  what  inputs  and  boundary  conditions 
must  be  recognized,  b)  what  outputs  are  expected,  and  c)  what 
COLISEUM  resource  functions  are  available.  The  ICD  and 
COLISEUM  resource  library  may  be  distributed  to  outside 
groups  so  that  COLISEUM  can  be  expanded  through  addition 
of  new  plasma  simulation  modules. 

The  COLISEUM  API  standardizes  the  definition  of  a  3-D 
plasma  problem  by  providing  strict  specifications  on  three 
categories  of  information:  surface  geometry,  surface 
properties,  and  sources.  The  manners  in  which  these  three 
categories  of  information  are  defined,  input,  stored,  and 
accessed  are  described  below. 

Surfaces 

Surfaces  are  modeled  in  finite-element  fashion  as  contiguous 
triangular  elements  joined  at  the  vertices  (nodes), 
COLISEUM  does  not  generate  3-D  geometries  or  surfaces; 
instead,  it  imports  them  from  other  software. 

Users  create  custom  geometries  using  almost  any  mainstream 
commercial  3-D  solid  modeling  package.  Then,  they  use 
finite  element  analysis  software  to  mesh  the  surface  of  their 
geometry  as  if  they  were  going  to  perform  a  structural 
analysis  using  thin  shells.  The  user  then  saves  the  meshed 
surface  file  in  ANSYS  format,  which  is  readable  by 
COLISEUM.  ANSYS  finite  element  format  was  chosen 
because  it  is  widely  supported  by  finite  element  packages. 

This  concept  of  separating  the  surface  geometry  definition 
from  the  plasma  calculation  has  proven  very  successful.  It 
greatly  reduced  development  time  and  cost  by  eliminating  the 
need  for  a  separate  surface  definition  module.  It  allows  users 
to  choose  which  software  to  use  in  defining  geometries.  And, 
users  can  import  into  COLISEUM  geometries  that  have 
already  been  defined  for  other  reasons  (structural,  thermal, 
etc.). 
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Surface  Properties 

The  user  may  provide  three  databases  in  conjunction  with  a 
surface  geometry:  a  component  database,  a  material  database, 
and  a  material  interaction  database. 

The  component  database  associates  specific  surfaces  with 
component  names  and  material  names.  These  associations  are 
established  by  using  a  component  number  which  is  specified 
in  the  ANSYS  file  using  integer  values  in  the  elastic  modulus 
field.  For  example,  the  component  database  may  specify 
component  number  34  as  component  name  “north_solar_cell” 
and  material  “quartz.” 

The  material  database  associates  component  names  with 
material  names  and  material  properties.  The  plasma 
simulations  RAY  and  PRESCRIBED_PLUME  require,  in 
addition  to  material  name,  molecular  weights,  and  charges  (in 
the  case  of  ions). 

The  material  interaction  database  contains  the  sputter  yield 
coefficients  and  sticking  coefficients  of  one  material 
interacting  with  the  other,  e.g.  between  Xe+  and  Kapton. 

Sources 

Sources  are  modeled  as  having  a  specific  velocity 
distribution,  fs(r,v,t) ,  that  is  a  function  of  position  on  the 
surface,  of  three-dimensional  velocity  space,  and  of  time: 

*»,  =  l\fs(r,v,t)d\ds  (1) 

sv 

Rather  than  specify  fs(rtv,t)  directly,  however,  three 

COLISEUM  resource  library  functions  are  provided  for  each 
source  type.  These  a)  give  the  distribution  of  velocities  at 
some  point  P  in  space  due  to  the  source,  b)  provide  a  random 
sample  from  the  source  velocity  distribution  at  the  surface,  or 


#  coliseum,  in 

# 

#  Load  the  Chamber  6  geometry, 

#  superinpose  the  2G0W  HET  plume 

#  calculation  by  SAIC,  calculate 

#  the  flux  and  net  sputtering  using 

#  ray  tracing,  and  save  the 

#  results  in  Tecplot . format . 

# 

c omp onen t_l o ad  component . txt 

material_load  xnaterial.txt  mat_mat.txt 

surface_load  ANSYS  Chambers .  ANS 

prescribed_plume_load  2DCIRC 
plume_SAIC_200W.dat  0,0  0.0  0.0  0.0  0.0  1.0 

ray  DEPOSIT  40 

surface_save  TECPLOT  Chamber6.dat 
FLTJXNORMAL .  XE+  SPUTTERRATE 


Fig.  2.  Sample  COLISEUM  command  file 


c)  update  the  source  to  be  valid  at  some  new  time,  t. 

This  method  is  extremely  descriptive  and  general.  Plasma 
simulation  modules  may  use  the  three  source  functions  to 
treat  the  source  distribution  function  in  various  ways.  For 
example,  using  the  first  function,  a  plasma  simulation  module 
can  be  written  to  treat  the  source  element  as  a  source  for 
geometric  ray  tracing.  Alternately,  particle  methods  can  use 
die  second  sample  from  the  velocity  distribution  and 
introduce  particles  randomly  over  the  full  element  surface. 
Therefore,  this  choice  of  source  definition  methods  gives 
COLISEUM  great  flexibility  by  enabling  a  wide  variety  of 
plasma  simulation  techniques  with  the  same  source  definition. 


Plasma  Simulation  Modules 

Two  plasma  simulation  modules  are  operational:  RAY,  a 
simple  ray  tracing  simulation  developed  by  AFRL,  and 
AQUILA,  an  unstructured  tetrahedral-mesh  PIC-DSMC 
simulation  developed  by  MIT.  In  addition,  Virginia 
Polytechnic  Institute  (VPI)  is  developing  a  PIC-DSMC 
simulation  which  uses  immersed  mesh  techniques.  The 
details  of  RAY  and  AQUILA  have  been  presented  before. 4>5 

This  paper  will  present  results  using  RAY  and  AQUILA. 


Execution  Sequence 

To  run  a  COLISEUM  case,  three  steps  are  typically  required: 
geometry  definition,  surface  meshing,  and  execution  of  the 
simulation. 

Geometry  Definition 

Users  generate  the  3-D  geometry  using  any  suitable  software. 
AFRL  uses  Solid  Works™.  Enclosed  geometries  such  as 
vacuum  chambers  are  typically  created  as  hollow  objects.  In 
these  cases,  it  is  convenient  for  the  geometry  definition 
software  to  have  cutaway-capability  so  that  the  interior  of  the 
enclosure  can  be  visualized. 

Surface  Meshing 

The  surface  meshing  process  may  also  be  performed  by  a 
variety  of  software  packages.  The  only  requirements  are  that 
the  software  must  produce  output  in  one  of  the  following 
formats:  ANSYS,  MGEN,  and  NASTRAN.  AFRL  is 
currently  using  COSMOS  DesignStar,  which  is  a  finite 
element  analysis  package  that  supports  ANSYS  output 
format. 

Execution  of  the  Simulation 

COLISEUM  executes  batch  commands  that  the  user  provides 
in  a  text  input  file.  The  commands  are  executed  sequentially 
as  they  appear  in  the  input  file.  Each  command  may  have 
some  number  of  parameters  separated  by  spaces  or  commas. 
A  sample  input  file  is  shown  in  Fig.  2. 
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Geometry  definition  typically  takes  approximately  6  hours  for 
medium-complexity  geometries.  Typical  run  times  for  low- 
fidelity  cases  (using  PRESCRIBED_PLUME  or  RAY)  take 
approximately  20  minutes  on  a  2  GHz  Intel  Pentium  4 
workstation.  Once  more  detailed  physics  are  incorporated 
with  plasma  algorithms  such  as  PIC-DSMC,  run  times  are 
expected  to  be  between  20  minutes  and  20  hours,  depending 
on  the  level  of  fidelity  and  on  the  initial  conditions. 

This  illustrates  a  key  feature  of  COLISEUM.  From  scratch,  a 
user  can  define  a  complete  three-dimensional  problem,  and 
generate  a  first  order  solution  all  in  less  than  one  workday. 
Then,  for  higher  fidelity  solutions,  the  problem  does  not  have 
to  be  redefined.  Since  the  plasma  simulation  modules  are 
interchangeable,  a  higher-fidelity  algorithm  may  be 
immediately  started  for  an  overnight  run. 

Acceleration  of  Particle  Simulations 

One  universal  difficulty  modeling  thrusters  in  vacuum 
chambers  is  the  long  computation  time  required  for  the 
system  to  relax  to  a  steady  state.  We  are  incorporating  several 
techniques  into  COLISEUM  to  mitigate  this  problem. 

There  are  three  main  phenomena  that  require  temporal 
resolution  during  the  plasma  simulation  of  the  thruster:  the 
particle  velocity,  the  electric  field,  and  the  collisions.  We 
must  resolve  the  particle  crossing  of  spatial  gradients. 
Although  these  simulations  are  relaxing  to  an  equilibrium 
state,  we  must  ensure  that  the  field  solution  does  not  suffer 
from  numerical  error  or  instability.  We  must  also  resolve  the 
collision  frequencies  of  the  various  particles.  In  this  section 
we  shall  quantify  these  time  scales  and  show  the  utility  of 
numerical  subcyeiing  for  some  particle  species. 

The  ion  flow  velocities  are  in  the  range  of  2-104  m/'s,  and  the 
neutral  thermal  velocity  is  near  200  m/s.  The  smallest 
simulation  volume  elements  are  0.01  m  on  a  side.  To  properly 
sample  spatial  gradients  the  particles  should  traverse  less  than 
one  third  of  an  element  during  a  particle  move.  Thus,  the 
particle  motion  time  limit  is  10‘7  for  ions  and  10'5  for  neutrals. 

In  the  PIC-DSMC  calculations  the  plasma  is  assumed  to  be 
quasi-neutral  with  electrons  subject  to  the  Boltzmann 
equation.  The  finite  difference  leap-frog  method  is  used  to 
advance  the  position  and  velocity  of  the  ion  particles.  We  can 
use  linear  theory  to  determine  the  stability  and  accuracy  of 
this  scheme  as  a  function  of  numerical  discretization  in  space 
and  time.  For  simplicity  consider  the  fluid  equations  for  the 
plasma 
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Assume  a  small  perturbation  to  the  equilibrium 
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Keeping  only  linear  terms  and  discretizing  with  finite 
differences  in  space,  j  Ax,  and  time  lAt,  we  find 
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We  can  combine  these  equations  to  eliminate  u  and  <b.  This 
results  in  the  dispersion  relation 

z-2  +  l/z  =  (/? /4)(w— 2  +  1/  w)-a(l-l/z) 
where  z  =  exp(-  imAt),  w  =  exp(zA'Ax),  a  =  vAt , 

f  A*\2 

.  For  the  highest  frequency 
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disturbance  the  wavelength  is  2Ax.  So  kmax  =  it/Ax,  and  w  =  * 
1.  The  result  is 

z2  -  (2  -  a  -  p)z  +  (l  -  a)  =  0 

The  solution  for  z  is  given  by 


Z  = 


1- 


a  +  p 


+  . 


2  J 


For  small  a  and  p,  the  magnitude  and  phase  are  given  by 


a 


and 


©  =  arctan 


(l-a) 

'-(] 

L-l 

{a  +  fi)/2) 

i2 

1/2  ’ 

1- 

(a  +  fi)/2 

Consider  Xe  with  Te  =  2  eV  and  Ax  =  1  cm.  The  magnitude 
and  phase  are  plotted  versus  the  temporal  discretization  in 
Fig.  3.  The  algorithm  remains  stable  for  At  less  than  5  usee. 
The  phase  error  becomes  significant  for  At  greater  than  0.5 
usee. 
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Fig.  3.  Magnitude  and  phase  of  the  shortest  wavelength 
numerical  mode 


The  following  formulas  were  used  for  the  collision  cross 
sections.  For  momentum  transfer  between  neutrals6 


2.117-10“ 


For  momentum  transfer  between  ions  and  neutrals6 


The  cross  section  for  charge  exchange  between  neutrals  and 
doubly  charged  ions  is6 

a  =  (-2.7038 ln(vre,)+35.006l)2  -10 -20m2. 


The  mean  time  between  collisions,  t,  is  given  by  the  inverse 
of  n<j<vtci>.  The  maximum,  neutral  density  in  these 
simulations  is  approximately  1018/m3,  Table  1  lists  the 
collision  times  given  the  typical  relative  velocities  between 
species,  A  simulation  time  step  of  less  than  10  usee  is 
necessary  to  properly  resolve  collisions. 


Table  1:  Collision  properties 


Collision 


Xe-Xe++  CEX 


6,8e-19 


8.2807-10- 


The  preceeding  analysis  indicates  a  situation  where  the 
phenomena  to  be  modeled  operate  on  two  disparate  time 
scales.  The  field  solution  and  ion  motion  require  a  time  step 
near  0,1  usee  while  collisions  and  neutral  motion  require  a 
time  step  near  10  usee.  For  computational  speed  we  have 
chosen  to  subeyele  the  ion  motion  and  field  solve.  Savings  in 
collision  computations  and  the  neutral  particle  push  lead  to 
significant  improvements  in  computation  time, 

A  flow  chart  of  the  subcycle  algorithm  is  shown  in  Fig,  4,  The 
main  loop  encompasses  the  procedures  to  advance  all  the 
particles  one  full  time  step  about  10  usee  for  these 
simulations.  Within  each  main  loop  the  fast  particles, 
typically  the  ions,  are  subeycled.  For  the  results  shown  below 
there  are  100  subcycle  steps. 

The  subcycle  loop  includes  the  time  advance  of  the  particle 
velocities  and  positions.  Fast  particles  are  injected  from 
sources.  The  E-field  is  also  updated  each  iteration.  Since  the 
code  presently  assumes  quasi-neutrality  with  Boltzmann 
electrons,  the  potential  is  obtained  directly  from  the  ion 
density.  Since  this  procedure  is  computationally  fast,  it  does 
not  significantly  impact  the  speed  of  the  subcycle  iterations. 
In  the  future  when  the  potential  is  calculated  from  Poisson’s 
equation,  the  solution  of  the  elliptic  equation  on  the  mesh  will 
modify  the  results  presented  here. 

After  completion  of  subcycling,  positions  and  velocities  of  the 
slow  particles  are  advanced  one  full  time  step.  Slow  particles 
are  injected  from  surface  sources.  Finally,  all  particles 
participate  in  DSMC  collisions. 


The  cross  section  for  charge  exchange  between  neutrals  and 
singly  charged  ions  is7 


cr  =  (-  23 .3  log10  (vrei )  + 142.21)  •  1 . 1 872  •  lO-20  m2 
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Fig*  4.  Flow  diagram  for  the  time  advance  with  subcycling. 

Results  and  Discussion 


Timing  Study 

A  set  of  simulations  were  run  to  quantify  the  effect  of 
subcycling  on  the  computation  time.  Neutrals  and  singly 
charged  ions  were  loaded  into  a  cube  with  1  m  sides.  The 
neutrals  were  given  a  density  of  1017  /m3  and  a  temperature  of 
300  K.  The  ions  were  given  a  density  of  1014/m3  and  a 
temperature  of  1 1600  K.  The  weights  were  set  so  there  were 
100,000  neutral  macroparticles  and  from  10,000  to  100,000 
ion  macroparticles.  The  time  step  for  a  normal  simulation  was 
0.1  usee.  The  time  step  for  a  subcycling  simulation  was  10 
usee  with  a  subcycle  time  step  of  0.1  usee.  Calculations  were 
run  with  and  without  collisions.  The  total  time  simulated  was 
1  msec. 

Fig.  5  shows  the  effect  on  computation  time  as  the  number  of 
ions  is  increased.  With  subcycling  most  of  the  time  is  spent 
moving  the  ions.  The  difference  with  and  without  collisions  is 
only  a  few  percent  causing  the  two  curves  to  overlap.  The 
zero  intercept  represents  the  time  required  to  move  the 
neutrals.  Without  subcycling  all  particles  are  moved  and 
collisions  occur  every  0.1  usee  time  step.  Comparing  the 
normal  curve  to  the  normal  curve  without  collisions, 
approximately  60%  of  the  time  is  taken  up  in  computing  the 
collisions.  The  result  is  that  even  as  the  number  of  ions  equals 
or  exceeds  the  number  of  neutrals,  the  subcycle  simulation 
still  runs  3  times  fasters  than  the  normal  simulation.  For  lower 
ion  particle  numbers  the  subcycle  algorithm  can  achieve  a 
speed  increase  of  15  or  more. 


Fig.  5.  Change  In  computation  time  as  a  function  of  the 
number  of  ions  In  the  timing  simulation. 

Vacuum  Chamber  Simulation 

To  illustrate  the  full  range  of  capabilities  of  COLISEUM,  a 
Hall  thruster  cluster8,9  was  simulated  using  both  RAY  and 
AQUILA.  A  cluster  of  four  thrusters  was  arranged  in  a  2  by  2 
matrix  in  AFRL  chamber  6.  The  thrusters  were  200  watt 
Busek  BHT-200-X3.  These  operate  with  discharge  voltage 
and  current  of  250  V  and  0.80  A  respectively.  The  anode  and 
cathode  mass  flows  are  8.5  seem  and  1.0  seem.  The  simulated 
chamber  and  thruster  configurations  are  shown  in  Fig.  6. 

RAY  calculates  fluxes,  sputter  rates,  and  deposition  rates  on 
surfaces.  These  results  are  shown  in  Fig.  7  and  Fig.  8. 
Unfortunately,  no  experimental  sputter  or  deposition  data  was 
collected  for  this  geometry. 

AQUILA  was  used  with  subcycling  to  calculate  local 
densities,  velocity  distributions,  and  plasma  potentials  on  the 
surfaces  and  in  the  plume  region.  The  density  of  ions  and 
neutrals  in  the  chamber  are  shown  in  Fig.  9  and  Fig.  10.  Fig. 
11  compares  the  calculated  current  density  at  0.5m  to 
experimental  measurements  made  at  MIT  and  AFRL. 


Conclusions 

Although  still  in  an  early  stage  of  development,  COLISEUM 
now  can  help  predict  ion  flux  and  equilibrium  net  sputtering 
and  deposition  rate  of  surface  materials  both  onboard 
spacecraft  and  in  laboratory  test  facilities.  COLISEUMS 
modular  architecture  is  allowing  rapid  expansion  of  its 
capabilities,  and  giving  users  flexibility  to  design  their  own 
geometries  and  choose  their  preferential  plasma  simulation 
method. 

Additional  work  for  the  future  includes  further  investigation 
of  the  re-sputtering  process,  further  validation  against 
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experimental  data,  and  construction  of  new  plasma  simulation 
modules  that  can  self-consistently  compute  plasma  expansion 
and  interaction  with  surfaces. 
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Fig.  11.  Charge  density  versus  angle  at  50  cm  from  the 
exit  plane  of  one  thruster. 


